Life and death of a cardiac calcium spark 



Michael D. Stern/ Eduardo Rios,^ and Victor A. Maltsev^ 

'Laboratory of Cardiovascular Science, Intramural Research Program, National Institute on Aging, National Institutes of Health, 
Baltimore, MD 21224 

^Department of Molecular Biophysics and Physiology, Rush University, Chicago, IL 60612 

Calcium sparks in cardiac myocytes are brief, localized calcium releases from the sarcoplasmic reticulum (SR) be- 
lieved to be caused by locally regenerative calcium-induced calcium release (CICR) via couplons, clusters of ryano- 
dine receptors (RyRs) . How such regeneration is terminated is uncertain. We performed numerical simulations of 
an idealized stochastic model of spark production, assuming a RyR gating scheme with only two states (open and 
closed). Local depletion of calcium in the SRwas inevitable during a spark, and this cotild terminate sparks by in- 
terrupting CICR, with or without assumed modtilation of RyR gating by SR lumenal calcium. Spark termination by 
local SR depletion was not robust: under some conditions, sparks could be greatly and variably prolonged, termi- 
nating by stochastic attrition-a phenomenon we dub "spark metastability." Spark fluorescence rise dme was not a 
good surrogate for the duration of calcium release. Using a highly simplified, deterministic model of the dynamics 
of a couplon, we show that spark metastability depends on the kinetic relationship of RyR gating and junctional SR 
refilling rates. The conditions for spark metastability resemble those produced by known mutations of RyR2 and 
CASQ2 that cause life-threatening triggered arrhythmias, and spark metastability may be mitigated by altering the 
kinetics of the RyR in a manner similar to the effects of drugs known to prevent those arrhythmias. The model was 
unable to explain the distributions of spark amplitudes and rise times seen in chemically skinned cat atrial myo- 
cytes, suggesting that such sparks may be more complex events involving heterogeneity of couplons or local propa- 
gation among sub-clusters of RyRs. 



INTRODUCTION 

Resting, mammalian cardiac muscle produces local cal- 
cium release events known as calcium sparks (Cheng 
et al., 1993). It is widely accepted that these events are 
caused by coordinated CICR (Fabiato, 1983) among 
local clusters ("couplons"; Stern et al., 1997, 1999) of cal- 
cium release channels (RyRs) , usually located in a dyadic 
cleft between the SR and the sarcolemma. This creates a 
paradox. If spark release is caused by local, regenerative 
positive feedback of CICR, what causes or allows it to ter- 
minate? Early on, it was believed that a time- and Ca^*- 
dependent inactivation of RyRs was involved (Stern and 
Cheng, 2004) , but studies of isolated RyRs of the cardiac 
isoform in lipid bilayers have not demonstrated a suffi- 
ciently powerful inactivation process (although it re- 
mains possible that such inactivation exists in vivo or is 
mediated by allosteric interactions, as suggested by Stern 
et al., 1999) . It is now generally accepted that local deple- 
tion of calcium in the terminal cisternae of the junctional 
SR (JSR) is involved in spark termination. A common hy- 
pothesis (Lukyanenko et al., 1998; Terentyev et al., 2002) 
is that when [Ca'*]jsR falls below a critical threshold (often 
suggested to be ^^50% of the resting level), sparks ter- 
minate because of deactivation of RyRs, either by direct 
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action of a lumenal calcium-sensing site on the RyR or 
mediated by the calcium-buffering protein calsequestrin 
(Casq) located in the JSR. 

In this paper, we examine this hypothesis semi- 
quantitatively by means of numerical modeling. Specifi- 
cally, we generate statistical distributions of spark 
properties from a minimal model of stochastic CICR 
in spatially extended RyR clusters, and try to determine 
whether, and under what conditions, observed spark 
statistics could be explained by such a model. We find 
that spark termination can be produced by local JSR 
depletion, with or without the hypothesized regulation 
of RyR gating by lumenal calcium. However, the robust- 
ness of spark termination depends on kinetic relation- 
ships among rates of RyR gating and of SR calcium 
recycling. In other words, the ability of depletion to 
terminate sparks is not absolute but depends on the 
kinetics of RyR gating and JSR refilling. An important 
consequence of this result is prediction of a novel regimen 
of sparks, dubbed "metastable sparks," which persist for 
a prolonged period before terminating by stochastic at- 
trition. Anything that makes the JSR refill more quickly 
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or causes a greater sensitivity of the RyR to calcium, or 
a slower RyR closure rate (longer open dwell time) , will 
favor the metastable spark regimen. 

MATERIALS AND METHODS 

Calcium sparks were simulated by an extension of a stochastic 
algorithm described previously (Stern et al., 1999). In brief, the 
algorithm generated a statistical ensemble of 10,000 realizations 
of a model couplon consisting in a square array of RyRs located 
on a 30-nm lattice in a dyadic cleft. Sparks were initiated by ran- 
domly opening one RyR in the couplon. Dyadic cleft [Ca^*] was 
calculated by a reaction-diffusion equation discretized on a two- 
dimensional 10-nm grid as described previously (Stern et al., 
1999). In addition, the differential system now inckided a JSR 
compartment, with instantaneous bviffering by Casq, and, option- 
ally, a chain of free SR (FSR) compartments without buffering. 
The JSR compartment (Fig. 1) was coupled to an infinite pool of 
FSR calcium through either a single-diffusion resistance (Fig. 1 A) 
or a chain of 100 compartments simulating the local FSR of 
the sarcomere (Fig. 1 B). The joint gating of the RyRs in a cou- 
plon was simulated by a Monte Carlo algorithm that produces 
an exact realization of the variable-rate, continuous-time Markov 
process driven by the evolving cleft calcium distribution (Stern 
etal., 1997). 

For this study, we confined ourselves to a minimal gating 
scheme in which a single RyR has only two states, closed and 
open, with no inactivation. The transition rates were given as a 
function of local cleft and lumenal calcium concentrations by: 

openingrate = k„[Ca^^^'2,ji{sloCf^ + [Ca^^^j^j^) 



closingrate - k^^, (1.2) 

where hr^ is the cooperativity of RyR opening by calcium, k„ is the 
opening rate constant of the RyR (dimension mM"*""* ms"'), and 
the parameter sloco determines the relative change in opening 
rate as lumenal calcium is decreased from resting. The usual value 



of sloco, 0.1 mM, implies an 11-fold decrease of opening rate as 
SR calcium falls from 1 mM to zero. For the simulations pre- 
sented here, hry, the cooperativity of RyR activation by cytosolic 
calcium, was taken to be 2, in keeping with many studies of iso- 
lated RyRs in lipid bilayer (see Qin et al., 2009, for a graphical 
summaiy of many results). We did not include allosteric cou- 
pling ("coupled gating"; Marx et al., 2001) among adjacent RyRs, 
as demonstrations of this phenomenon in RyR2 have been few 
and controversial. 

To compare simulations with experimental spark fluorescence 
data, the total RyR flux from each couplon was taken as source 
for a spherically symmetrical diffusion and buffering computa- 
tion, and convolved with an empirical confocal point-spread 
function, as described previously (Santiago et al., 2010). Histo- 
grams of the distribution of observable in-focus spark properties 
(e.g., rise time, amplitude) were constructed from the ensemble 
of detected sparks out of 10,000 individual spark simulations gen- 
erated in parallel. Sparks were considered to be "detected" if they 
satisfied the criteria that peak F/Fq > 1.2, with a peak rate of rise 
>0.2 ms"'. 

To examine further the mechanisms of spark initiation and ter- 
mination, we computed the behavior of a simple, deterministic 
two-compartment continuum model of a couplon ("toy model"). 
Further details of the model structures, parameters, and imple- 
mentation are given in the Appendix. 

Online supplemental material 

The online supplemental material includes a description of struc- 
ture and implementation of the stochastic model, RyR gating es- 
timation of Dup, and deterministic simulations of sparks/blinks. 
Fig. SI shows simulated central spark fluorescence using a variety 
of different endogenous buffer parameters taken from published 
papers, assuming the same, decaying release flux through a fixed 
nvimber of open RyRs. The calcvilated spark amplitude varies 
widely depending on the assumptions. Fig. S2 shows JSR calcium 
computed from the full, deterministic simulation using a con- 
stant square-wave release flux, compared with a single-compart- 
ment model using the value of D^p estimated from steady-state 
diffusion as described in the supplemental text. Fig. S3 shows that 
peak fluorescence is reached early, regardless of the length of 
RyR opening. Even for large rates of JSR refilling, depletion 

Figure 1 . Structure of the model. The model 
consists of 10,000 independent couplons, each 
with a dense square cluster of RyRs in a dy- 
adic cleft, with an associated JSR release ter- 
minal. The JSR is refilled by intra-lumenal 
diffusion from calcium in the FSR (SR uptake 
compartment), which is assumed to be an in- 
finite compartment at resting calcium CaSRo 
because isolated sparks are being considered. 
Two models of the refilling connection were 
considered. In the simplest (A) , the JSR refills 
through a single diffusional resistance, char- 
acterized by the refilling rate constant D^p. As 
discussed in the Results, it is difficult to rec- 
oncile the observed refilling rate with simple 
diffusion theory. We therefore constructed a 
second version of the model (B) in which the 
connection is made through a tube of local 
FSR, whose length and diameter are chosen to 
match the observed steady-state diffusion re.sis- 
tance (characterized by time constant Tmi) and 
the observed volume fraction of FSR. For the 
standard parameters and a true half-sarcomere 
length of 1 pm, the effective SR tube length 
is 1.99.5 pni. 
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causes fluorescence to decline. Table SI lists typical model param- 
eters. Videos 1 and 2 show the RyR gating dynamics of a normal 
spark and a metastable spark, respectively. The online supple- 
mental material is available at http://www.jgp.org/cgi/content/ 
full/jgp.201311034/DCl. 



RESULTS 



How big is a spark? 

To determine whether it is feasible to compare simu- 
lated spark amplitudes with experimental ones, we per- 
formed a series of deterministic spark fluorescence 
computations, each using the same calcium source, but 
cissuming four different systems of cytosolic buffers, each 
of which had been used in a published spark modeling or 
analysis paper (Smith et al., 1998; Soeller and Cannell, 
2002; Santiago etal., 2010; Shkryl et al., 2012). As shown 
in Fig. SI, the peak amplitude of the simulated sparks 
ranged over an order of magnitude depending on which 
buffering system was assumed. We conclude that no ab- 
solute comparison can be made between the amplitudes 
of simulated and observed sparks unless the parameters 
of cytosolic diffusion and buffering have been mea- 
sured in the same preparation from which the sparks 
were recorded. Henceforth, we concentrate on statisti- 
cal properties of spark populations that are indepen- 
dent of absolute amplitude. 

Estimation of JSR refilling 

Modeling sparks requires estimation of the rate at which 
calcium from the FSR network refills the JSR release ter- 
minal. The detailed architecture and topology of the 
cardiac FSR are not well known. The simplest model as- 
sumes that the JSR refills from an infinite FSR pool 
through a single effective diffusion resistance. This re- 
sistance-scaled to JSR volume-may be characterized by 
the parameter D^p (dimension ms"^), which is the rate 
constant of JSR refilling that would pertain in the ab- 
sence of Casq. Estimation of Aip from first principles 
turns out to be problematic. 

A simple model is to treat the FSR as a fine-grained 
random network occupying a volume fraction fvfi,- 
{'^3.5%; Shannon et al., 2004) of the cell. By a theorem 
of stereology, this is also the fraction of any surface 
crossed by the FSR lumen. As a result, the diffusion 
equation for FSR calcium is formally the same as that 
in the cytosol, except that boundary fluxes are reduced 
by the factor fafi,- If we assume spherical geometry, with 
the JSR occupying a sphere of radius rjsrand volume Vp„ 
we can estimate D,,p as the steady-state diffusional con- 
ductance (flux over concentration difference) between 
rjs, and infinity, divided by the JSR volume. Expressing 
r^„ in terms of the fractional volume /i7,5rOfJSR and the 
length and radius of the sarcomere gives (see supple- 
mental text) : 



3I/3 42/3 j^^^^j)^ 

Uup = 



fVjs. 



2/3; 



2/3 „ 4/3' 

isarc 



(1.3) 



where D^^ is the calcium diffusion coefficient in the SR 
lumen, estimated to be 0.06 |im^/ms (Picht et al., 2011). 

If we assume that the JSR occupies 0.35% of the vol- 
ume of a cylindrical sarcomere of length 2 pm and di- 
ameter 0.8 pm, this gives D^p = 0.70765 ms ^ Although 
this calculation was based on a steady-state diffusion 
gradient. Fig. S2 shows that a single-compartment 
JSR simulation using this value, including 30-mM Casq- 
binding sites, reproduces well the JSR depletion and 
refilling of a full, deterministic spark/blink simulation. 
The problem is that, in both simulations, the refilling 
phase (^10-ms time constant in the presence of Casq) 
is still much too rapid when compared with experi- 
mentally observed blinks, which have time constants 
of ^200 ms in rabbit (Picht et al., 2011), or experimen- 
tal spark restitution, which recovers approximately as 
a 90-ms exponential in rat (Sobie et al., 2005; Ramay 
et al., 2011). Alternatively, we may estimate Aip from 
these observed time constants, by assuming that the ob- 
served time constant t is the exponential time constant 
of a one-compartment JSR with Casq, linearized about 
the resting state (see supplemental text) : 



Dup = 



{kdcsq^ + 2casrO kdcsq + CasrO^ "jz 



(1.4) 



This gives values for Z)up in the range of 0.05-0.15 ms \ 
considerably less than that estimated above from the 
random network model. 

A more detailed simulation of FSR diffusion has been 
given by Picht et al. (2011), in which the FSR network 
was modeled as linear segments connecting JSR termi- 
nals in a three-dimensional grid. We ran their model 
with an assumed 20-ms square-wave release current. Al- 
though the resulting JSR calcium (not depicted) had a 
nonexponential time course, the time to e-fold recovery 
was ^30 ms in the presence of Casq, corresponding to 
A,/, = 0.3 ms \ Picht et al. (2011) hypothesized that the 
much slower observed blink recovery was caused by 
continuing RyR release on the tail of a spark. As we 
show below, this interpretation is problematic. 

In summary, there is a discrepancy between the ob- 
served rates of blink refilling and spark restitution versus 
simple models of intra-lumenal diffusion in the SR. 
Presumably, the discrepancy will be accounted for ei- 
ther by a more accurate measurement of the intra-SR 
calcium diffusion coefficient or by further details of the 
nano-anatomy of the SR, for example, a constriction 
or bottleneck in the connection of the JSR to the FSR 
(Brochet et al., 2005) . In the meantime, we will assume 
that Aip lies somewhere in the range of 0.05-1 and 
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consider the consequences of these possible vahies for 
spark termination. 

Lumenal depletion is unavoidable 

Most experimental studies have estimated that only 
^50% of JSR calcium is released during a spark. In the 
Appendix we calculate how much lumenal [Ca^^] would 
be reduced, in steady state, by an indefinitely prolonged 
opening of Wopen RyRs. Fig. 2 shows that the steady-state 
JSR depletion produced by even a single, open RyR 
would be >50% even at the highest Aip, corresponding 
to a refilling time of '^10 ms, much faster than observed. 
Obviously, the large number of RyRs that participate 
in a spark would be expected to produce much more 
severe depletion, especially with a more realistic refill- 
ing rate. Therefore, major JSR depletion is inevitable 
unless a spark were terminated very early by a very pow- 
erful inactivation process. This suggests that the degree 
of local JSR depletion ("blink") during a spark may have 
been underestimated by confocal microscopy. 

The meaning of "spark duration" 

It is often stated that the rise time of spark fluorescence 
at the center of the spark is an approximation of the 
duration of calcium release during the spark. This is not 
true for the rapidly declining release flux from a deplet- 
ing JSR terminal. As shown in the supplemental text. 
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Figure 2. Steady-state calcium concentration in the JSR that 
would result from holding open different numbers of RyRs, based 
on Eq. A1.18, as a function of the JSR refilling rate constant D^p. 
The computation assumes that a single RyR has a unitary current 
of 0.35 pA (Meji'a-Alvarez et al., 1999) at a resting [Ca^'^JjsR of 1 
mM, and takes into account the series resistance of the dyadic 
cleft. The parameter D^p is the reciprocal of the time constant 
for JSR refilling by diffusion from the FSR that would pertain in 
the absence of Casq. The relationship of Aip to refilling time in 
the presence of Casq is given in Eq. 1.4. For comparison, Dup = 1 
would give a refilling time of ^10 ms, much faster than observed. 
The observed rate of spark restitution in rat (Sobie et al., 2005; 
Ramay et al., 2011) would correspond to D^p of ~0.1, and blink 
recovery time in rabbit (Picht et al., 2011) would be even lower 
The implication is that opening of even a single RyR is capable of 
depleting JSR calcium by >50% under any conceivable assump- 
tion about intra-lumenal diffusion. Note that this steady-state cal- 
culation does not depend on the presence or any assumptions 
about Casq. 



rapid depletion of the JSR causes the amplitude and 
timing of the fluorescence peak to be controlled by the 
amount of calcium in the JSR, independent of the dura- 
tion of the release. From a mechanistic point of view, it 
would be better to define spark duration as the time 
elapsed between the opening of the first RyR that trig- 
gers the spark and the closing of the last RyR causing 
extinction of CICR. Fig. 3 A shows the time courses of 
the number of open RyRs; JSR free calcium and dye 
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Figure 3. Lack of relationship between spark fluorescence rise 
time and extinction time. (A) Time course of fluorescence F/Fq, 
the number of RyRs, and JSR free Ca^* for a representative cal- 
cium spark. (B) Histograms (1-ms bin width) of fluorescence rise 
time and spark extinction time. (C) Scatter plot of fluorescence 
rise time versus spark extinction time. 
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fluorescence for a representative, stochastic, and simu- 
lated spark; as well as histograms of fluorescence rise 
time and spark extinction time from "detected" sparks 
out of 10,000 stochastic sparks attempted. There is clearly 
no similarity between these two measures of spark dura- 
tion, and, as shown by the scatter plot (Fig. 3 B) , there 
is very little relationship between them. Unfortunately, 
spark extinction time is not reliably observable. Once 
JSR depletion has occurred, the tail of calcium release 
flux is likely to be lost in the noise inherent in inverting 
the spark fluorescence signal. 

Implicit in the histogram of Fig. 3 is that sparks do 
self-extinguish in a reasonable time. In other words, 
under the assumed conditions, JSR depletion provides 
a sufficient spark termination mechanism without re- 
quiring any RyR inactivation process. As shown in Fig. 4, 
this termination does not require the dependence of 
RyR gating on lumenal [Ca^^]. Spark extinction time 
was only modestly increased when lumenal activation 
was turned off. This shows that the interruption of posi- 
tive feedback of CICR among RyRs within a couplon-by 
reduction in RyR flux as the JSR becomes locally de- 
pleted-is a sufficient mechanism to terminate sparks (as 
also shown by two models recently published: Gillespie 
and Fill, 2013, and Laver etal., 2013). 

The robustness of spark termination, however, is 
strongly affected by the rate of refilling of the JSR 
from the neighboring FSR. The sparks in Figs. 3 and 4 
were generated under conditions when D^p was set to 
0.15 ms^^, as inferred from empirical restitution and 
blink relaxation rates by Eq. 1.4. If Z^^p is increased to 
0.7 ms^^-the value predicted for a random network 
model of the SR-then a substantial fraction of sparks 
fails to terminate early, and continues to "smolder" for 
a prolonged period, forming a straggling tail on the 
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Figure 4. Histograms of spark extinction time with and with- 
out the lumenal activation factor multiplying the RyR opening 
rate (Eq. 1.1). The lumenal activation factor provides about a 
10-fold variation in opening rate over the range of JSR calcium 
concentrations during a spark, but its absence causes only a mod- 
est prolongation of RyR opening, and sparks still terminate based 
entirely on the direct effect of JSR depletion to reduce RyR flux 
and CICR loop gain. 



histogram of extinction times (Fig. 5) . We refer to this 
predicted behavior as "spark metastability." 

An obvious (but insufficient) explanation for the fail- 
ure of these sparks to terminate would be that continu- 
ous, rapid refilling of the JSR prevents [Ca'^JjsR from 
being reduced to the "threshold" for spark termination. 
As Fig. 2 shows, even a modest number of open RyRs 
can severely deplete JSR calcium, even at the highest 
refilling rate; as shown below, this occurs rapidly on 
a time-scale short compared with spark duration. An- 
other explanation is needed. 

A toy model of the couplon 

To examine the spark extinction mechanism more 
closely, we constructed a deterministic "toy" model of a 
single couplon. We treat the couplon as a "bag" of RyRs 
sharing a single, common cleft compartment whose 
calcium is in instantaneous steady state. The number of 
open RyRs at a given instant is treated as a continuous 
variable. The flux through an open RyR is proportional 
to the difference between [Ca^*]jsR and [Ca^^Jcieft- 
[Ca"^] cleft is determined by the balance between aggre- 
gate RyR flux into the cleft and diffusion out of the 
edges of the cleft: we estimate [Ca^^Jdeft as the average 
value of [Ca'^] in a circular, thin (two-dimensional) 
cleft if the total RyR flux were deposited at the center. 
[Ca'*] cleft then depends only on the flux, the calcium 
diffusion coefficient, and the height (but not the diam- 
eter) of the cleft (see Appendix). The instantaneous 
"number" of open RyRs varies at a rate equal to the dif- 
ference between the opening rate and closing rate as 
determined by the kinetics of the RyR gating scheme, 
Eqs. 1.1 and 1.2. 

With these simple assumptions, the couplon becomes 
a deterministic dynamical system with two state vari- 
ables, n-open and ca^r, governed by a pair of differential 
equations that appear rather daunting once all the sim- 
ple pieces above have been assembled (see Appendix): 

(1.5) 

dca. 
dt 

(fc, + m.) I ^ ' 



SnccL^FhD — Snca^FhD + L 



- + 2ca,.D.,Fv,.. - 2ca,D, Fv,. 



(1.6) 



dt 



I 8 n ca. ca,. — 8nca, )FhD + ca, L 
{?>nca., - ?>nc.CL\FhD + Lru. 



We consider first Eq. 1.6 for ?i„^^„, with ca^r held con- 
stant. Taking the RyR calcium cooperativity hry to be 2, 



Stern et al. 261 



1500-1 



(A 

Q. 
CO 



0) 

E 

Z 



1000- 



500- 



,Dup=0.15 




0 100 200 

Spark Extinction Time (ms) 



300 



a "bifurcation point" at time 0) by spontaneous activa- 
tion by background calcium, it will rapidly, regenera- 
tively rise until most RyRs are open. Once this happens, 
cUsr will fall, raising the CICR threshold until it rises 
above w„^(,„. At this point, CICR is no longer regenera- 
tive, and n^pen will fall because of independent, sponta- 
neous closure of RyRs. This will allow the depleted 
JSR to begin refilling, so the threshold will fall again, 
following n„j,f„ downward. If the falling threshold re- 
mains above the live value of n„pg„, n„pen will converge to 
the stable, closed fixed point and remain there until a 
spontaneous opening triggers another spark. 

The condition just described corresponds to the 
green curves in Fig. 7, which shows the coupled dynam- 
ics of the deterministic model Eqs. 1.5 and 1.6 {n„pe„ in A, 
casr in B, and the phase-plane plot of n^pen vs. ca^r in C) . 
Those curves were computed with D^p = 0. 1 , a value that 
corresponds to the restitution time of 90 ms observed 
by Ramay et al. (2011). Suppose, however, we assume 
Aip = 0.7, the value derived from modeling the FSR as a 
fine, random network. In this case, our simulations re- 
vealed a rather complex behavior illustrated by the red 
curves in Fig. 7. What has happened is that the more 
rapid refilling of the JSR allows the falling CICR thresh- 
old to catch up with the falling Wo^„„ so that CICR again 
becomes regenerative. The system can "latch up," be- 
coming permanently trapped in state with a small num- 
ber of RyRs open, with their release flux balanced by 
steady refilling of the JSR, and the resulting cleft cal- 
cium sufficient to maintain their (average) open state. 
As shown in the phase-plane plot (Fig. 7 C), this indi- 
cates the creation of a new, stable fixed point of the 
combined system, onto which the trajectory converges 
instead of returning to the fully closed state. This cor- 
responds, in the toy model, to the "metastable spark" 
observed in the stochastic model. As discussed further 
in the Appendix (especially Fig. Al), the presence of 
this new fixed point is a steady-state phenomenon, de- 
pendent on the relationship between the steady-state 
values of n^pin and the steady-state depletion of the 
is raised above this threshold (shown in Fig. 6 A as JSR, which depends on but not on Casq. However, 

Figure 6. Threshold behavior 
of the toy couplon at con- 
stant SR calcium. (A) Dynam- 
ics of the "toy model" of the 
g couplon, when the variable 

ca„ ([Ca^*]jsR) is held at a 
fixed value. There is a sharp 
threshold value of M„^e„ above 
10H which CICR will regener- 

ate until almost all RyRs are 
open, and below which it 
will fall exponentially to near 
zero. (B) The CICR thresh- 
old is a steep function of JSR 
calcium, allowing CICR to 
cease regeneration when the 
JSR is depleted. 



Figure 5. Spark extinction times (5-ms bin width) assuming two 
values of the JSR refilling rate parameter D^^, with the RyR gat- 
ing rate constants held fixed at k„ = 100 mM"' ms"' and A„,„ = 0.1. 
Increasing D^^ to the value predicted from modeling the FSR as 
a diffuse, random network causes some sparks to "smolder" for a 
prolonged period rather than extinguishing promptly. 



this equation has three fixed points (i.e., equilibrium 
values of w„^j„ at which its rate of change vanishes) . One 
stable fixed point describes a resting condition with 
nopen = 0, a tiny value because of openings stimulated by 
resting background calcium. Another stable fixed point 
describes the condition where n^pen ~ n^y and most RyRs 
are open, kept open by the large cleft [Ca'^] caused by 
their currents. Between these two lies an unstable fixed 
point, describing a threshold (at a given casr) for regen- 
eration of CICR. As shown in Fig. 6 A, if n^g„ is started 
below the threshold (green curve) , it will fall to the 
closed fixed point, whereas if it is above (red), it will 
increase (initially exponentially) to converge to the 
fully open state. The value of the threshold depends in- 
versely on casr, as shown in Fig. 6 B. 

We can now describe the sequence of events when 
Eqs. 1.5 and 1.6 are coupled. Initially, the JSR is fully 
loaded, so the threshold for CICR is very low (Fig. 6 B) , 
reflecting the fact that RyRs are highly sensitive and 
RyR current is high, favoring regeneration of CICR. If 
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whether a metastable spark actually forms depends criti- 
cally on the relative rates of the various processes. Rela- 
tively small changes in parameters can affect whether 
the trajectory is captured by the new equilibrium point. 

In the more realistic, stochastic model, this trapped 
state is not permanently stable. Random fluctuations 
allow the system point to "leak" out of the attraction 
basin of the equilibrium point, permitting metastable 
sparks to extinguish slowly by a process of stochastic 
attrition, giving the long tail on the histogram of spark 
duration. Fig. 7 D shows phase plots of the ensemble 
average ri„^„, and [Ca^^J^jr of 10,000 stochastic sparks in 
stable and metastable regimes. In the metastable regi- 
men, sparks pause near a "latch point" before slowly re- 
turning to the resting state by stochastic attrition. This 
is a fundamental difference between deterministic and 
stochastic (more realistic) behavior. 

What the toy model shows is that spark metastability 
is a kinetic matter, depending on the quantitative rela- 
tionship among the rate of refilling of the JSR and the 
gating rates of the RyR. Anything that makes the JSR 
refill more quickly, the RyR more sensitive to calcium 
(including removal of regulation by lumenal calcium) , 
or the RyR closure rate slower (open dwell time longer) 
will favor the metastable spark regimen. This is a matter 
of practical concern, because the prolonged release 
from metastable sparks may permit CICR to propagate 
to adjacent couplons (not considered in this paper) 
causing diastolic calcium release events that may result 
in triggered arrhythmias. In fact, mutations of Ryr2 
and CASQ2 having these kinetic effects are known to 
cause catecholaminergic polymorphic ventricular tachy- 
cardia in humans and animal models (Jiang et al., 2004; 
Cerrone et al., 2005, 2007; Liu et al., 2006; Chopra et al., 
2007; Paavola et al., 2007). Likewise, the arrhythmias 



Figure 7. Free dynamics of the toy model for 
stable (green) and unstable (red) values of 
the JSR refilling rate. (A) "Number" of open 
RyRs as a function of time; (B) [Ca^"^]jsR ; 
(C) phase-plane parametric plot of ra„^„ versus 
ca^ (D) Phase-plane plot of number of open 
RyRs and [Ca^*]jsR in the ensemble average 
of 10,000 stochastic sparks generated by the 
full model for stable and unstable parameters. 
Green, stable: Tm = 6.5 ms (D„f = 0.15); red, un- 
stable: Tflu = 1 ms (D^p =1). Compare with C. 

that may accompany heart failure, and ischemic and 
reperfusion injury, may be explained by changes (phos- 
phorylation, oxidation, nitrosylation) that affect gating 
of the RyR in ways leading to spark metastability (Eisner 
et al., 2009). The conditions for spark metastability are 
also critically dependent on the diffusion coefficient of 
free calcium in the dyadic cleft, as discussed below. This 




Spark Extinction Time (ms) 

Figure 8. Attempt to "pharmacologically treat" spark metastabil- 
ity. The black curve shows a histogram (log ordinate) of spark 
extinction times, assuming RyR gating rates estimated from the 
lipid bilayer data of Guo et al. (2012) and the JSR refilling rate es- 
timated from spark restitution in rat (Sobie et al., 2005) . As shown 
by the tail extending to long durations, these parameters give un- 
stable sparks. The green curve shows the effect of doubling both 
the opening and closing rate constants of the RyR, thereby leaving 
the "calcium sensitivity" defined by p^ unchanged. This substan- 
tially shortens the unstable tail, because it accelerates RyR gating 
relative to the rate of refilling of the JSR. The red curve shows 
the additional effect of doubling the concentration of Casq, slow- 
ing down JSR depletion and refilling by increasing buffering. The 
peak of the distribution moves to the right, as expected, but the 
tail is actually shortened because JSR refilling is slower relative to 
RyR gating, reducing spark metastability. 
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Figure 9. An overlay of "^50 sparks in the metastable regimen 
defined by the gating parameters from Guo et al. (2012) and the 
JSR refilling rate inferred from observed spark restitution. The 
metastability is only revealed by the presence of low amplitude 
embers after the peak of the spark. 

variable is not generally amenable to modification, but 
its value in vivo is poorly known. 

"Treating" spark metastability 

Treatments that redress the kinetic balance should be 
helpful for such arrhythmias. Fig. 8 shows spark extinc- 
tion time histograms (log ordinate scale) for a model 
with RyR gating rates taken from a two-state approxima- 
tion to the lipid bilayer data of Guo et al. (2012) (see 
supplemental text for a discussion of the issues in mak- 
ing this estimate) and assuming = 0.15. As seen from 
the black curve, these parameters (taken from observa- 
tions in different kinds of experiments) actually give 



sparks that are somewhat metastable. The green curve 
shows the effect of simultaneously doubling the RyR 
opening and closing rate constants, keeping the po as a 
function of [Ca^^Jcieft unchanged. This kinetic change 
substantially reduces the tail of metastable sparks. The 
red curve shows the additional effect of doubling the 
concentration of Casq, thereby slowing JSR depletion 
and refilling. Increasing Casq delays the peak of the 
histogram as expected, but the tail actually moves left- 
ward as the refilling rate decreases compared with the 
RyR closure rate, reducing spark metastability. This of- 
fers a possible explanation for the observation that drugs, 
such as carvedilol and flecainide, that increase the RyR 
closure rate (decrease mean open time) have been ef- 
fective in preventing calcium waves and catecholamin- 
ergic polymorphic ventricular tachycardia in animal 
models (Watanabe et al., 2009; Hilliard et al., 2010; Zhou 
et al., 2011). 

The paradox of restitution 

As demonstrated by Picht et al. (2011) and confirmed 
by our simple spherically symmetric spark/blink compu- 
tations (see above and supplemental text) , the observed 
rates of blink recovery (^^200 ms in rabbit ventricle) 
and spark restitution ('^90 ms in rat using repeating 
sparks generated by micro-dose ryanodine; Sobie et al., 
2005; Ramay et al., 2011) are not easily accounted for by 
simple models of calcium diffusion in the FSR lumen. 
Possible explanations are the presence of tortuosities or 
constrictions in the connection of the JSR to the FSR, or 
an error in the estimation of the intra-lumenal calcium 
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Figure 10. Experimental spark amplitude dis- 
tribution in cat atrium cannot be accounted for 
by model. (A) Amplitude distribution of in-focus 
sparks in skinned cat atrium, replotted from 
Shkryl et al. (2012). (B) Distributions of peak 
F/Fo of simulated sparks with several different 
couplonsizes. (C) Amplitude (F/Fq-I) distribution 
of sparks simulated for a 7 x 7 couplon assuming 
different values of the RyR opening rate con- 
stant ko. For the highest value, sparks are unstable, 
as shown by the extinction time distribution in 
the inset. 
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diffusion rate. Picht et al. (2011), however, interpreted 
this discrepancy as evidence that there is a continuing 
tail of calcium release during the recovery from a blink, 
and that the time constant of recovery is limited by the 
decay of this release rather than by the rate of refilling 
from the FSR. Our analysis of spark metastability above 
makes this explanation problematic. If [Ca^^JsR recovers 
while some RyRs are still open, CICR will not extin- 
guish: the spark would be metastable. On the assump- 
tion that metastable sparks are not a normal condition, 
we take this as evidence against Picht's interpretation 
and favor the assumption that greater limitations on 
intra-SR diffusion are actually present. However, in the 
interest of open-mindedness, we must at least consider 
the hypothesis that "normal" spark/blink pairs are, in 
fact, metastable. For sparks that are terminated by a JSR 
depletion signal, the time-to-peak of spark fluorescence 
depends almost entirely on JSR volume, the number of 
RyRs recruited, and the resistance to efflux from the 
cleft. A metastable spark would only be recognized by 
a trailing "ember" (Gonzalez et al., 2000) because of 
continuing release through open RyRs, fed by calcium 
returning from the FSR. Fig. 9 shows an overlay of 
50 metastable sparks (F/Fq at the spark center). For 
these simulations, refilling rate was low but metastability 
was produced by raising the calcium sensitivity of the 
RyR. The resulting embers, although clearly visible in 
the noise-free simulation, might be difficult to detect 
under the conditions of confocal recording. Note also that 
the embers tend to be oscillatory. So called "sub-sparks" 
have in fact been noted after major sparks (Santiago 
et al., 2013) but have usually been attributed to recruit- 
ment of sub-clusters of RyRs. In our simulations, couplons 
are single crystalline arrays of RyRs, without satellite 
groupings, so oscillatory embers are a purely dynamical 
phenomenon (see Videos 1 and 2). 

Spark amplitude distribution 

In chemically skinned cat atrial cells, Shkryl et al. 
(2012) measured the distribution of amplitudes and 



rise times for a population of in-focus sparks, using a 
rapid two-dimensional confocal scanner combined 
with piezo z movement to locate the center of the spark 
in the z direction. Their data (replotted in Fig. 10 A) 
show that amplitudes are broadly distributed, with a 
long tail in the direction of higher amplitudes. For 
comparison. Fig. 10 B shows simulated spark ampli- 
tude distributions for several couplon sizes. For these 
simulations, we made the assumption that JSR volume 
is proportional to the number of RyRs in the couplon. 
That is, we assumed that the cell contains 1.4 x lO*" 
RyRs (Bers and Stiffel, 1993) and that the whole-cell 
JSR occupies 0.35% of cell volume, partitioned equally 
among couplons. The larger couplons therefore pro- 
duce larger sparks mainly because they have access to 
larger stores of JSR calcium. However, for any given 
couplon size, the amplitude distribution is rather nar- 
row, with a sharp upper limit constrained by JSR con- 
tents. Fig. 10 C shows that for fixed couplon size, the 
amplitude distribution also depends on the RyR open- 
ing rate constant ft„, which determines how many and 
how rapidly RyRs will be recruited and, therefore, the 
speed with which the JSR contents will be released. 
Even if RyR sensitivity is reduced to the point that 
sparks become highly stochastic, the amplitude distri- 
bution preserves the sharp upper cutoff. We conclude 
that the long tail on the spark amplitude distribution 
observed by Shkryl et al. (2012) cannot be accounted 
for by a homogeneous distribution of couplons with 
spark termination caused by local depletion. To ac- 
count for the data would require either heterogeneity 
of couplon sizes or the presence of more complicated 
events involving local propagation of CICR between 
couplons. This is further confirmed by the joint distri- 
bution of amplitudes and rise times (Fig. 11). The ob- 
servations show a positive correlation between rise time 
and amplitude. For simulated sparks, the relationship 
is negative for the fundamental reason that a longer 
rise time implies a slower release of approximately the 
same amount of JSR calcium. 
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Figure 1 1 . The experimental relation- 
ship between spark rise time and am- 
plitude is not explained by the model. 
(A) Two-dimensional joint histogram of 
rise time and amplitude of sparks from 
skinned cat atrium, taken from Shkryl 
etal. (2012). (B) Analogous two-dimen- 
sional scatter plot of simulated sparks 
from the stochastic model, which fails 
to reproduce the positive correlation 
seen in the experimental data. 
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Residual JSR calcium 

Fig. 12 A shows the distribution of the nadir of JSR free 
calcium during a spark, with or without lumenal modu- 
lation of RyR gating. The fraction of calcium remaining 
in the JSR after a spark is quite low, as has been found 
in other models (Sobie et al., 2002; Ramay et al., 2011; 
Laver et al., 2013) . This is at variance with blink observa- 
tions showing only ^50% depletion, as well as with the 
common idea that JSR calcium falls to a "threshold" of 
^50% at which point lumenal deactivation terminates 
the spark. 

Although the limited resolution of confocal micros- 
copy probably explains the relatively shallow blinks 
(Kong et al., 2013) , we attempted to determine whether 
different conditions might allow spark termination with 
a higher true JSR calcium nadir. Reducing the RyR 
open time to 1 ms to allow RyRs to close rapidly once 
the CICR threshold is reached, and reducing the RyR 
calcium sensitivity, can increase the average residual cal- 
cium. However, the nadir of JSR calcium did not ap- 
proach 50% unless RyR sensitivity was so low that sparks 
became very stochastic, with low amplitude and a very 
broad multimodal distribution of residual JSR calcium 



rather than a "threshold" (Fig. 12 B). In this condition, 
only 7% of spark attempts ignited a detectable spark. 
The difficulty in reproducing the hypothesized lumenal 
threshold behavior is caused by the fact that a substan- 
tial number of open RyRs will deplete the JSR very rap- 
idly on a time scale comparable to RyR gating. This is 
a fundamental constraint caused by the assumed ratio 
of RyRs to JSR volume and the magnitude of the RyR 
unitary current. We varied the parameters of cytosolic 
and lumenal calcium regulation of RyR gating in an 
attempt to determine whether there was a regimen 
in which lumenal "threshold" behavior was observed. 
The only manipulation that was able to recreate such a 
threshold was to make the RyR closure rate (i.e., open 
dwell time) strongly and very nonlinearly dependent on 
lumenal calcium, effectively "slamming shut" the RyRs 
when lumenal calcium falls to the threshold (Fig. 12 C). 
As there is no experimental evidence for such strong 
lumenal regulation of RyR closure rate, our simula- 
tions support the belief that the apparent 50% residual 
calcium in "blinks" is caused by limitations of confocal 
imaging of the JSR terminal, which is at the resolution 
limit of light microscopy. The issue remains unresolved. 
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Figure 12. The nadir of JSR calcium 
during a simulated spark is low unless 
RyR closing rate is strongly regulated 
by lumenal calcium. (A) Distributions 
of the nadir of [Ca^*]jsR of simulated 
stochastic sparks with (green) or with- 
out (red) the lumenal activation fac- 
tor in Eq. 1.1. Lumenal deactivation 
when the JSR is depleted contributes to 
termination of the spark at somewhat 
higher [Ca^*]jsR, but the nadir is quite 
low in both cases, which is not com- 
patible with the idea that the spark is 
shut off when lumenal calcium falls to 
a "threshold" of ~50%. (B) Attempt to 
reproduce the 50% nadir by lowering 
the RyR mean open time so that sparks 
will extinguish promptly once CICR is 
no longer regenerative, and reducing 
the opening rate constant. This re- 
quired RyR sensitivity so low that only 
7% of attempted sparks were detect- 
able, and they were very stochastic and 
either abortive or low amplitude. (C) A 
well-defined 40% nadir was produced 
by making the RyR mean open time 
strongly and nonlinearly dependent on 
[Ca^*]jsR, thereby forcing RyRs to close 
below the threshold. 
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Figure 1 3. RyR activation is incom- 
plete during a spark. (A) Histogram 
of the number of RyRs opened dur- 
ing a spark for various RyR calcium 
sensitivities. Even at the highest 
level, when sparks become meta- 
stable, not all RyRs open. (B) En- 
semble averages of the number of 
open RyRs and JSR calcium during 
the spark. The JSR calcium pool 
depletes rapidly during the period 
when RyRs are still being recruited 
by local propagation of CICR in the 
couplon, causing the spark to termi- 



Time (ms) nate before full recruitment. 



however, as there is strong experimental evidence 
(namely, agreement between depletion measured in 
blinks and in more global images) suggesting that 
the limiting depletion is well measured locally (Zima 
etal., 2008). 

Extent of RyR recruitment 

Somewhat surprisingly, for the larger couplon sizes not 
all RyRs open even during a robust spark. Fig. 13 A 
shows histograms of the number of activated RyRs for 
several RyR opening rate constants, as in Fig. 12 C. The 
explanation is that lumenal calcium falls at a rate compa- 
rable to the rate at which activation of RyRs propagates 
across the couplon (Fig. 13 B), so that CICR regenera- 
tion may fail while many RyRs have not yet been 
opened. The propagation of RyR activation is shown in 
Videos 1 and 2. 

Diffusion of free calcium in the dyadic cleft 

The resistance to exit of calcium from the dyadic cleft 
depends on the local diffusion coefficient of free cal- 
cium (Dca), which is not direcdy known. In free solu- 
tion, Dca depends on temperature, having a value of 
0.748 |iVms at 23°C and 1.05 at 37°C (Hrabetova et al., 
2009) . This value is reduced by viscosity and "tortuosity" 
in cytosol; values of 0.25-0.35 have been assumed in 
spark modeling studies (Smith et al., 1998; Soeller and 
Cannell, 2002; Laver et al., 2013). The dyadic cleft is 
largely occupied by RyR protein, so it is likely that 
further tortuosity effects exist. Tadross et al. (2012) 
found that the calcium nano-domain in the vicinity of 
an L-type channel pore is 10 times stronger than pre- 
dicted, implying a highly restricted diffusion environ- 
ment (i.e., Dca ~ 0.07). For most of our simulations, we 
used Dca = 0.15. 



Fig. 14 shows the distributions of spark extinction 
times assuming RyR kinetics derived from the results of 
Quo et al. (2012) as described above, and using values 
of Dca of 0.7, 0.35, 0.15, or 0.07 pVms. It is apparent 
that the degree of spark metastability is markedly in- 
creased by restricted calcium diffusion in the cleft; for 
the value of Dca implied by Tadross et al. (2012) , 75% of 
sparks fail to terminate within 1 s. This problem cannot 
be dealt with by rescaling the model: cleft geometry 
is constrained by ultrastructure, JSR contents are con- 
strained by ultrastructure and biochemical measure- 
ments of CSQ, and RyR kinetics and conductance are 
constrained by lipid bilayer studies, so CICR loop gain is 
determined once Dca has been specified. This places a 
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Figure 14. Spark extinction time distributions as a function of 
assumed free calcium diffusion coefficient in the dyadic cleft. In 
these simulations, couplons consisted of 7 x 7 RyRs, with k„„ = 
0.117 ms"', Ca50 = 10 [iM, JSR volume of 0.35% of cell volume 
divided among 20,000 couplons, and refilling time constant of 
90 ms in the presence of 30 mM Casq. Only detected sparks that 
terminated within 300 ms are included in the histograms. 
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premium on understanding the extent to which diffu- 
sion and channel properties are modified by the in situ 
environment of the cleft. 

DISCUSSION 

Calcium sparks in cardiac muscle are generally under- 
stood to be caused by locally regenerative CICR within a 
cluster of RyRs. How such regeneration, once ignited, is 
terminated is a matter of uncertainty and controversy. 
A commonly espoused idea is that local depletion of 
calcium in the JSR causes deactivation of RyRs-by a lu- 
minal calcium-sensing site on either the RyR itself or 
Casq-once JSR calcium falls below a critical threshold 
(Lukyanenko et al., 2001; Zima et al., 2008). We per- 
formed numerical experiments on a simplified model 
of stochastic CICR to test whether this mechanism can 
terminate sparks in the absence of any inactivation of 
RyRs. Using a minimal two-state RyR gating scheme, the 
initiation of sparks by a single RyR opening and their 
termination in response to JSR depletion could be re- 
produced. The model was deliberately simplistic to study 
the mechanisms in a reductionist manner. Therefore, 
we made no attempt to fit quantitatively the experimen- 
tal data. To do so would be unwise, in any case, because 
crucial parameters of RyR gating, diffusion within the 
SR, and buffering in the cytosol show wide discrepan- 
cies among studies. 

The most important finding in our study is that the 
ability of depletion to terminate sparks is not absolute 
but depends on the kinetics of RyR gating and JSR re- 
filling. When RyR calcium sensitivity is high, RyR open 
dwell time is long, and/or intra-lumenal diffusion of 
calcium in the SR is rapid or buffering in the JSR is 
weak, sparks can fail to terminate cleanly. Instead, they 
enter a period of prolonged, stuttering, often oscilla- 
tory, low grade release. Such release is ultimately termi- 
nated randomly by stochastic attrition, leading to a long 
tail on the distribution of spark durations. We refer to 
this phenomenon as "spark metastability." 

Such prolonged RyR opening has the potential to 
trigger calcium release from adjacent couplons, result- 
ing in propagated diastolic CICR or waves, which can 
trigger life-threatening arrhythmias. Spark metastability 
might be difficult to recognize when it is produced by 
high RyR sensitivity with normal refilling of the JSR, 
as the trailing calcium release from the depleted JSR 
might be inconspicuous. Zima et al. (2008) demon- 
strated extremely prolonged sparks when RyRs were 
partially blocked by tetracaine. They interpreted this as 
evidence that refilling of the JSR was rapid enough to 
prevent [Ca^^JjsR from being depleted to the "thresh- 
old" at which release would be terminated. The obser- 
vation is similar to our spark metastability, except that 
release flux appeared to remain constant for hundreds 
of milliseconds, as though only a single channel were 



open or RyRs gated synchronously. This phenomenon 
appears to be specific to tetracaine, as no similar pro- 
longed openings were seen when RyR unitary current 
was reduced by Tris (Guo et al., 2012), which should 
have similar effects if the tetracaine effect is caused by a 
reduced efflux of JSR calcium through partially blocked 
RyRs, as proposed by Zima et al. (2008). Long-lasting 
openings of the kind seen by Zima might well be the 
result of a more complicated gating scheme than the 
minimalist one we deliberately restricted ourselves to. 
It is possible that tetracaine produces long-lasting open 
states or subconductance states (as ryanodine does) 
that have not been studied, that there is some form of 
inactivation that appears in vivo although not seen in 
bilayer and, most intriguingly, that cooperativity caused 
by aJlosteric interactions might result in collective modes. 
These possibilities can be explored in future work 
using the same model structure but open up a vast ter- 
ritory of parameters and gating schemes that are not 
considered here. 

Somewhat unexpectedly, the ability of JSR depletion 
to terminate sparks did not require the modulation 
of RyR gating by lumenal calcium. In the absence of 
such modulation sparks still terminated. This indicates 
that the interruption of positive CICR feedback by de- 
creased RyR flux is a sufficient and probably dominant 
mechanism of spark termination. A similar conclusion 
has been reached, indirectly, from the effects of reduc- 
ing RyR unitary currentwith Tris (Guo et al., 2012), and 
from recent modeling studies by Laver et al. (2013) 
using a model of a couplon with slighdy more realistic 
geometry than ours. Even when lumenal modulation of 
RyR opening rate was present, we found litde evidence 
for the existence of a deterministic "threshold" of SR 
calcium at which sparks would terminate. Unless RyR 
closure rate was made strongly dependent on lumenal 
calcium, sparks had a very low nadir of JSR calcium. 
Similarly strong depletion of JSR calcium has been 
found in other models (Sobie et al., 2002; Ramay et al., 
2011; Laver et al., 2013). This suggests that the shallow 
depletion during "blinks" is an artifact of blurring by 
confocal microscopy. This has been demonstrated re- 
cently by Kong et al. (2013), who used a deterministic 
spark diffusion model (very similar to the one we used 
in Fig. S3) to fit observed spark/blink pairs, without as- 
sumptions about the underlying release mechanism. 
They found that the optimal fit required marked JSR 
depletion (masked by blurring), early decay of release 
flux, and RyR opening lasting well past the peak of fluo- 
rescence, as produced in our stochastic simulations. 

The role of JSR depletion in terminating sparks, 
and the associated phenomenon of spark metastability, 
places a premium on understanding the rate of refilling 
of the JSR. Diffusion in the SR lumen has been mea- 
sured by Picht et al. (2011) using fluorescence photo- 
bleaching recovery, obtaining a value of 0.06 pm^/ms. 
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which is an order of magnitude less than free calcium 
diffusion in aqueous solution, and approximately one 
fifth of the values usually assumed to pertain in cytosol. 
Their model, which assumed that there is a concentra- 
tion of 27 mM of calcium-binding sites on Casq in the 
JSR, estimated JSR refilling rates at least three times 
faster than observed blink recovery rates. Our naive 
computation, assuming a random FSR network extend- 
ing directly to the JSR, gave refilling rates several times 
faster yet. Picht et al. (201 1) rationalized this by assum- 
ing that there is continuing release during the blink re- 
covery phase. As we have shown, this can only be true 
if sparks are metastable. Alternatively, one must con- 
sider the presence of an additional diffusion resistance 
between the FSR and JSR, or a much higher concentra- 
tion of buffer in the JSR. Brochet et al. (2005) estimated 
in rabbit ventricular myocytes that there is an average 
of 4.5 discrete, narrow connections attaching the JSR 
release terminal to the FSR network. The concentration 
of Casq in the JSR is also poorly known, with esti- 
mates ranging from 10 to 120 mM of Ca-binding sites 
(Murphy et al., 2011). We therefore assumed the pres- 
ence of a phenomenological diffusion resistance con- 
necting the JSR with the distant FSR pool, characterized 
by a refilling rate parameter D^p. In the simplest form of 
the model, we treated this as a single resistance between 
the JSR and the infinite pool of SR calcium far from the 
spark. Because Picht et al. (2011) showed in their model 
that SR depletion during a blink can extend for several 
microns, we also constructed a version of the model 
in which the local FSR is represented as a long tube 
(Fig.l B) that is able to supply some of its calcium to the 
JSR during the course of a single spark. This resulted in 
somewhat larger spark amplitudes as expected (not de- 
picted) but did not change the qualitative behavior of 
the model. 

As discussed in the supplemental text, we attempted 
to estimate the parameters of our two-state RyR gating 
scheme from lipid bilayer measurements made in "cell- 
like salt solutions" (Guo et al., 2012). These resulted in 
sparks that were somewhat metastable, so most of our 
simulations were done using lower values of the open- 
ing rate constant k„ (equivalently, higher EC50 of the 
open probability curve). Laver et al. (2013) made sto- 
chastic simulations of sparks based on measured pa- 
rameters and obtained sparks that terminated promptly. 
This results from two factors. For unknown reasons, 
their lipid bilayer data in rat showed an EC50 an order 
of magnitude higher than that found by Guo et al. 
(2012) . Second, although most of the geometry of their 
sarcomere was based on observations, they assumed 
that the JSR filled from the FSR through a single 30-nm 
"nexus," which has an effect similar to our imposition of 
an empirically constrained low value of D,i/<- 

By studying the analogue of spark metastability in a 
simple, deterministic "toy" model of the couplon, we 



established that spark metastability is fundamentally a 
kinetic phenomenon, determined by the relationship 
between the refilling rate of the JSR, which controls 
spark refractoriness, and the gating rates of the RyR, 
which determine how quickly the open probability of 
the RyR diminishes in response to JSR depletion, pre- 
venting re-ignition. Because uncontrolled CICR can have 
life-threatening consequences, we considered whether 
it might be possible to "treat" spark metastability by alter- 
ing this kinetic balance. To demonstrate the kinetic na- 
ture of the effect, we increased both the opening and 
closing rate constants of the RyR by the same factor, 
leaving the po as a function of calcium unchanged. This 
gready reduced spark metastability (Fig. 8) , which was fur- 
ther improved by slowing JSR refilling by increasing 
Casq from 30 to 60 mM. Zhou et al. (2011) studied 14 (3 
blockers used to treat ventricular arrhythmias in heart 
failure and found that only one, carvedilol, could sup- 
press calcium waves. They showed that it acted by de- 
creasing the open dwell time of the RyR, which is 
equivalent to increasing the closing rate fto,„. ft had no 
effect on the opening rate constant (which might nor- 
mally be thought of as the "calcium sensitivity" of the 
RyR) and actually increased the total number of open- 
ings. However, because the closing rate controls the 
rate at which RyR openings respond to JSR depletion, 
as well as the rate of stochastic attrition, increasing the 
closing rate has the most effect on spark metastability, 
explaining why carvedilol is most effective in treating 
calcium-overload arrhythmiais (Poole-Wilson et al., 2003). 
Flecainide, another antiarrhythmic drug that reduces 
RyR open time, has also been shown to be uniquely pro- 
tective against calcium waves and triggered arrhythmias 
(Watanabe et al., 2009; Hilliard et al., 2010). 

An analysis similar to our toy model was presented by 
Hinch (2004), who examined the asymptotic properties 
of a simplified Markov model of sparks, leading to quasi- 
deterministic analytic results parallel to what we have 
done in the "toy." Hinch examined a somewhat differ- 
ent regimen (basically that of Sobie's "sticky cluster" 
model; Sobie et al., 2002), which favors bistable behav- 
ior (flat top sparks) with a narrow stochastic transition 
between open and closed. The phenomenon that we 
call "spark metastability" was demonstrated, using a 
null-cline analysis similar to our Fig. Al (Appendix), 
specifically as an effect of removing alios teric coupling, 
although the generality of the effect as a property of 
relative kinetic rates was not explored. 

A model is actually most informative when it is unable 
to fit the data. Figs. 10 and 1 1 show that in at least one 
model system-chemically skinned cat atrial myocytes- 
the distribution of in-focus spark amplitudes hais features 
that cannot be explained by our model. The fact that 
the observed amplitude distribution has a long tail that 
is positively correlated with spark rise time cannot be 
accounted for by a model in which spark termination is 
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caused mainly byJSR depletion in a homogeneous set 
of couplons. One can speculate that this might be evi- 
dence of heterogeneity of couplon properties, or, more 
radically, that events counted as "sparks" actually involve 
local propagation to more than one couplon. This lat- 
ter possibility seems especially plausible because the 
atrium, lacking t-tubules, depends on radial propagation 
of CICR to excite the center of the fiber. In service of this 
propagation, there may be small, nonjunctional RyR 
clusters with spacings as close as 100-150 nm (Franzini- 
Armstrong, C, personal communication) that might 
appear as a single "spark" source. This invocation of het- 
erogeneity may also be appropriate to account for an- 
other observation by Shkryl et al. (2012) that is not readily 
consistent with any of the termination mechanisms dis- 
cussed here, namely, that sparks of different rise time are 
essentially superimposable (in fluorescence and its rate of 
rise) up to their peak time. 

Summary and future work 

By means of numerical simulations, we have demon- 
strated that it is possible to achieve termination of cal- 
cium sparks in a model with RyRs having only two states, 
with no inactivation or adaptation process. Termination 
results from local depletion of JSR calcium, which is an 
inevitable consequence of the high conductance of the 
RyR channel and the architecture of the SR. Depletion 
terminates sparks by interrupting the positive feedback 
of CICR by reducing calcium flux into the dyadic cleft, 
and also by reducing the open probability of the RyR by 
removal of activation by lumenal calcium, although the 
latter mechanism is not essential. However, termination 
of sparks by these mechanisms is not absolute and can 
be greatly prolonged if the kinetic relationship among 
gating rates of the RyR and refilling of the JSR is un- 
favorable, a phenomenon we dub "spark metastability." 
This has the potential to contribute to the development 
of life-threatening triggered arrhythmias. Spark meta- 
stability is critically dependent on the calcium-depen- 
dent opening and closing rates of the RyR, the resistance 
to diffusion within the SR lumen, the concentration of 
Casq in the JSR, and the diffusion coefficient of free 
calcium in the dyadic cleft, none of which is known with 
precision. In particular, more evidence is needed re- 
garding restriction of diffusion within the dyadic cleft, 
which our simulations (Fig. 14) show would have a pro- 
found effect. 

Using the best estimates of these factors, sparks from 
couplons with '^BO RyRs may, in fact, be metastable. 
This indicates either that metastable sparks do occur in 
vivo, or that further experimental work is needed to bet- 
ter determine the conditions at the couplon. Particular 
uncertainty attaches to the rate of refilling of the JSR by 
diffusion within the SR lumen, as observed rates are sev- 
eralfold slower than might be expected from simple 
models. Restrictions to diffusion of calcium in the cleft 



need further study, as recent estimates are incompatible 
with spark termination in large couplons. 

The simulations shown here assumed a cooperativity 
of 2 for calcium activation of the RyR, as that is usually 
found in lipid bilayer in "cell-like" solutions with high 
lumenal-side [Ca'^]. However, it is difficult to explain 
the low resting spark rate under these conditions. Our 
qualitative findings are not altered by using higher co- 
operativity for RyR opening (not depicted) , but the pre- 
cise conditions for spark metastability may be different. 
One possible mechanism to explain high cooperativity 
in vivo would be the existence of allosteric interactions 
among RyRs in dense clusters (Stern et al., 1999). We 
have not included allosteric interactions in spark simu- 
lations because evidence for its existence is sparse, but 
it needs to be studied in future simulations, together 
with the fact that couplons have been found to consist 
of smaller sub-clusters (Baddeley et al., 2009). 

Finally, we must mention the most important idealiza- 
tion in our model: the fact that sparks are treated as 
isolated and independent. It is conventional to say 
that sparks are the elementary events of excitation- 
contraction coupling, from which whole-cell calcium 
transients are formed by superposition. In fact, how- 
ever, calcium release from one couplon will influence 
its neighbors in at least three ways: (1) cytosolic calcium 
will diffuse to activate, or at least bias, RyRs; (2) SR cal- 
cium will be drawn away from other couplons by dif- 
fusion within the lumen; and (3) SR calcium will be 
increased by diffusion through the cytosol followed by 
uptake by SERCA. Each of these effects has different 
spatio-temporal characteristics. The concept of CICR 
metastability needs to be extended to the whole cell, 
which is a hierarchical array of RyRs, a kind of couplon 
of couplons. Such extension goes beyond the scope of 
this paper. 

APPENDIX 

Construction of the "toy" model 

We first estimate the average free calcium in the dyadic 
cleft when a calcium current i is entering the cleft. We 
assume that the distribution of calcium in the cleft is in 
instantaneous steady state. The cleft is treated as a cylin- 
der of radius Rand (small) height h and analyzed in two 
dimensions. At the edge of the cleft, a boundary value 
of caO, equal to bulk cytosolic [Ca^*] , is assumed. This is 
a crude approximation but good enough for our quali- 
tative purposes. The calcium current i entering the cleft 
via RyRs is taken to all be deposited at the center. 

In this approximation, free calcium as a function of 
radius will satisfy the two-dimensional Laplace equation 
in polar coordinates, with circular symmetry, except at 
the origin, which is singular because of the delta-func- 
tion source: 
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d ca 
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dca 
dr 



(A1.7) 



The solution satisfying the boundary condition is easily 
verified to be: 



ca = a( log R - log r) + cao. 



(A1.8 



Differentiating and setting r = R, we find that the flux of 
calcium per unit area exiting at the edge of the cleft is 
given by: 

aD 



f- 



R 



(A1.9) 



from which the constant a can be determined by requiring 
that the total calcium flux out of the cleft be equal to the 
calcium entering via the current i This gives the solution: 



i{\ogR-logr) 

ca = — + cao, 

AnFhD 



(ALIO) 



where F'^s, the faraday. 

Integrating the above expression over the area of the 
cleft gives the average [Ca^*] in the cleft, which we will take 
cis the calcium seen by RyRs on their cytosolic face: 



calcium gradient across the pore, such that the unitary 
current under resting conditions {ca„ = ca„o, <ca> = cao) 
is iuo- 
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casro — cao 

Eliminating <ca> between these two equations gives: 

87t(casr — cao^Fh luO Tlopen D 



SncasroFhD - SncaoFhD + iuor, 



(A1.13) 



The total calcium concentration in the JSR, including in- 
stantaneous buffering by Casq, is: 



casrt = casrl — h 1 

kdaq + CasT 



(A1.14) 



casrt changes caused by calcium leaving through RyRs 
and being replenished from the distant SR pool through 
a single diffusion resistance: 



dcasrt , . i 

[casrO — CasrjDup 



dt 
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Performing the differentiation and solving for the rate 
of change of free calcium in the JSR gives: 



(ca)-- 



SnFhD 



- + cao. 



(Al.ll) 



Note that this expression depends on the height of the 
cleft but not on its radius. If we had assumed that the 
source current were uniformly distributed over the area of 
the cleft, we would have obtained a parabolic expression 
for cleft calcium, giving an average value differing from 
the above by a factor of two, which is immaterial for our 
qualitative purposes. Recall that, in the full, stochaistic 
model, the exact distribution of cleft calcium is calculated 
dynamically at the position of each RyR 

The current i is taken to enter through 'n„pg„ RyRs, each 
of which has linear permeation proportional to the free 
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dt 



Inserting the above-derived expression for i gives the 
ca^r dynamical Eq. 1.5 shown in the main text: 

(A1.17) 

dca. 
dt 
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Figure A1 . Phase-plane sta- 
bility analysis of stable (A) and 
unstable (B) sparks in the "toy 
model." The thin dashed black 
line plots the stationary bound- 
ary for TC„p„, and the thin ma- 
genta line plots that for ca^^ In 
A, the stationary curves do not 
intersect, so there is no global 
fixed point. Increasing _D„^ to 0.7 
(B) moves the magenta curve so 



that it intersects the black curve 
in two fixed points, stable and 
unstable. For the parameters of 

the red curve, the trajectory becomes trapped on the stable fixed point. Increasing Casq (blue curve) allows the trajectory to escape the 
fixed point without altering the location of the fixed point, showing that spark metastability is a matter of kinetics. 
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Setting the left-hand side to zero and solving for casr 
gives the steady-state value that JSR calcium would 
achieve if Wo^,,,, RyRs were kept open indefinitely, which 
is plotted in Fig. 2, for various values of D^f. 

(A1.18) 

^(Sncasro^ — Sn caocasro^DapFhvp^ + Ak cmhiuonopn^jD + casroDup Lo nopmVpr 
{ij!)n casro — Sncao)Dup Fhvpr + iKhiuo nop^^D + Dup iuo iup^Vpr 

We assume that the continuous variable n^pcn changes in 
the same way as the ensemble average for RyRs activated 
by a cytosolic-side calcium <ca>: 

- = (ca)*'' ke,(^nry - nopen)(^slOCO + CUsr) - komtlopen. (Al .19) 



dt 



where h,y is a Hill coefficient giving the cooperativity of 
RyR activation (taken to be 2 in the numerical examples 
in the text) . 

Inserting the above-derived expression for <ca> gives 
the dynamical Eq. 1.6 for n„pc„ shown in the main text: 



dViOpm 

dt 



(A1.20) 
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Phase-plane stability analysis of the toy model 

Here, we further clarify the steady-state and kinetic as- 
pects of spark metastability using the deterministic toy 
model. The model provides an approximation to spark 
metastability as trapping of the system trajectory on a 
stable fixed point, an aspect of the steady-state behavior 
of the system. Fig. Al shows how to locate such a point 
by a graphical analysis of the phase plane of versus 
ca^rfor stable (A) and metastable (B) parameter values 
(as in Fig. 7, but plotted with logarithmic ordinate for 
clarity) . In each panel, the dashed thin black line shows 
the curve (nullcline) along which dn„pc„/dt is zero. For 
points above the curve, will increase, whereas for 
points below, it will decrease. This curve is the same in 
both panels; its location does not depend either on Z)„^ 
or on the absolute scale of RyR gating rates or the pres- 
ence of Casq. The thin magenta line shows the nullcline 
along which dcasr/dt vanishes; i.e., it plots the steady- 
state depletion of the JSR. This is the same information 
as in Fig. 2 but plotted as a function of n„pg„ for fixed Z)„^. 
The location of this curve depends on D^^p but not on 
RyR gating rates or Casq. A fixed or stationary point of 
the full system will exist wherever these two curves inter- 
sect. For the parameters in Fig. Al A, where Z)„^ =0.1, 
there is no intersection, so no "latch point" exists and 
no spark metastability is possible. The system trajectory 



(Fig. Al A, green curve; same as in Fig. 7 C except for 
the logarithmic axis) returns exponentially to the rest- 
ing {n„pin = 0) state. In Fig. Al B, Ai/<= 0.7. This raises the 
magenta curve so that it now intersects the dashed Wo^^,, 
stationary curve in two places, creating a pair of fixed 
points, stable (bottom right) and unstable (top left). 
The de novo appearance of these equilibria as a func- 
tion of a parameter is an example of a saddle bifurca- 
tion. However, the existence of a stable fixed point does 
not guarantee spark metastability, which is determined 
by whether the trajectory is actually captured by the 
point. That is a matter of kinetics. For the particular 
parameters in Fig. Al B, the red trajectory is captured 
and spirals into the fixed point where it will remain. 
However, increasing Casq from 27 to 120 mM, which 
does not affect static stability (i.e., the location of the 
equilibrium point), allows the trajectory (blue curve) to 
circumvent the equilibrium point and return to the 
resting state. The same effect can be produced by in- 
creasing both the opening and closing rates of the RyR 
proportionately (not depicted). The possible role of 
Casq buffering to stabilize sparks was pointed out by 
MacLennan and Chen (2009) , but they interpreted it in 
terms of steady-state buffering preventing SR calcium 
from reaching a hypothetical threshold for store-over- 
load release. 
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